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Radiative shocks in galaxy formation. 

I: Cooling of a primordial plasma with no sources of 

heating. 
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ABSTRACT 

We use a 1-D Lagrangian code which follows both a gaseous and a dark component to 
study the radiative shocks that appear in the evolution of spherical scale- free perturba- 
tions in an Einstein-de Sitter Universe. The detailed behaviour of the shock depends 
on whether the radiative cooling is dominated by bremsstrahlung or line cooling. 
Bremsstrahlung is the main energy loss mechanism for systems with circular velocity 
V c > 100 km s _1 . In this case, we can reproduce the kinematics of the shock and of 
the cooling wave to a high degree of accuracy with a simple analytic model. When 
line cooling dominates, the shock front can be unstable to oscillations. The period 
and amplitude of the oscillations increase in time as the universe expands. For such 
systems, our analytic model provides only a rough estimate of the mean evolution. 
We believe that now we fully understand the effect radiative cooling has in the shocks 
that appear in spherical models for galaxy formation. 
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1 INTRODUCTION 

The origin of galaxies is still an open question in mod- 
ern cosmology. The standard picture of structure formation 
supposes that the nonlinear objects we observe today grew 
through gravitational insta bility from a spectrum of lin- 
ear primordial fluctuations (Peebles 1967). The cosmic mi- 
crowave background radiation (C MB) can severely constrai n 
such structure formation models ( |White, Scott fc Silk 1994 ) , 
but observations with high angular resolution are required 
to put direct constraints on scales approaching those rele- 
vant to galaxy formation. Cosmogonies dominated by cold 
dark matter (CDM) are the current standard, but there are 
few unambiguous tests of their predictions so far on galactic 
scales (Blumenthal et al. 1984; Davis et al. 1985; Ostriker 
1993; Ostriker and Steinhardt 1995). 

In a cold dark matter scenario structure forms hierarchi- 
cally, which means that small objects form first and merge 
to form bigger ones. As the gas falls into the potential well 
of the forming halo it passes through a shock in which its 
bulk kinetic energy is converted into thermal energy. In the 
case when the gas can cool efficiently it will flow further 
in, getting very dense. This cold and dense gas in the in- 
nermost parts of halos provides suitable raw material for 
star formation. The important role of radiative cooling in 
structure formation was first pointed out independently by 



Binney (1977), Rees & Ostriker (1977) and Silk (1977), and 



the way in which it leads to galaxy formation at the centres 
of dark matter halos was first described and modelled by 
White & Rees (1978). If we want to understand the forma- 
tion of galaxies it is essential to study how the gas cools and 
settles in the central regions of dark matter halos. 

Rapid improvements in computer technology have 
played an important role in the theory of structure forma- 
tion, since numerical simulations can provide detailed quan- 
titative predictions to be compared to real data. Until quite 
recently only the gravitational clustering of the collisionless 
component could be followed (e.g. Davis et al. 1985). Over 
the last few years, however, it has become feasible to intro- 
duce gas dynamics and cooling/heating processes in simu- 
lations of galaxy formation (e.g. Cen 1992; Katz, Hernquist 
& Weinberg 1992; Navarro & White 1993; Steinmetz 1996). 
Such three-dimensional simulations require a large amount 
of CPU-time, and as a result exhaustive parameter stud- 
ies are difficult. Moreover, numerical artifacts often compli- 
cate interpretation of the results, and a clear picture of the 
physics behind the evolution of such simulations is not easy 
to build. 

A very useful alternative is provided by semi-analytic 
models which combine an algorithm for predicting the abun- 
dance and merging of dark-matter halos (or alternatively a 
high resolution cosmological N-body simulation of this pro- 
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cess) with cooling arguments to prescribe galaxy condensa- 
tion at halo centres, with a simple star formation law, and 
with parametrised models for supernova feedback and for 
spectrophotometric evolution; such models are ideal for pa- 
rameter studies and for direct comparison with observation, 
and their quantitative predictions can agree reasonably well 
with a wide range of data (e.g. White & Frenk 1991; Kauff- 
mann, White & Guiderdoni 1993; Cole et al. 1994; Kauff- 
mann, Nusser & Steinmetz 1997). 

On the other hand, semi-analytical models are only as 
good as the simple models they use for the various phys- 
ical processes involved. One-dimensional numerical simula- 
tions are a potential tool to improve such models which do 
not suffer from many of the inconveniences of their three- 
dimensional counterparts. However, very little has been done 
in this direction since the pioneering work of Larson (1974). 
The formation and growth of plane-symmet ric pa ncakes has 
been studied by Shapiro & Struc k-Mar cell ( 1984 ) and more 
recently by Anninos & Norman (1994); Thoul & Weinberg 
studied dissipational spherical collapses onto local density 
maxima of a primordial G aussian density perturbation field 
( rhoul fc Weinberg 1995 ) and the effects of a ba ckground 
of ultr aviole t radiation (Thoul & Weinberg 1996); Haiman 
et al. (1996) adapted the code built by Thoul and Wein- 



berg to investigate radiative cooling by molecular hydrogen. 
The present paper lies within this general framework and 
focusses on the study of the radiative shocks that appear 
in the evolution of spherical scale-free perturbations in an 
Einstein-de Sitter universe. It is organized as follows: in Sec- 
tion 2 we present our basic model for formation of a galaxy, 
as well as the lD-Lagrangian code we use to study its evo- 
lution; in Section 3 we test the code against some existing 
solutions; in Section 4 we present the results of our simu- 
lations of galaxy formation; finally, in Section 5 we present 
physical explanations for the behaviour observed in the sim- 
ulations. 



2 MODEL AND NUMERICAL METHOD 

The observed rotation curves of spiral galaxies and the or- 
bital velocities of their satellites seem to indicate that the 
surrounding dark halos have an approximately isothermal 
profile p ~ r~ 2 (e.g. Casertano & Van Gorkom 1991; Zarit- 
sky & White 1994). An appropriate 1-D model for halo for- 
mation should build up such a halo thro ugh i nfall of ma- 
terial from an expanding universe. Gott (1975) and Gunn 
(1977) were the first to construct a model of this type. They 
considered accretion onto a central seed from an otherwise 
unperturbed universe, and showed that the resulting density 
profile had p oc r~ 9 ' 4 for an Einstein-de Sitter background 
cosmology. By allowing a wider varie ty of scale-free initial 
perturbations Fillmore & Goldreich (1984) showed that a 
range of power law expone nts ( including —2) could be ob- 
tained. White & Zaritsky ( |l992| ) (thereafter WZ) extended 

They showed 



these self-similar infall models in two ways, 
how to embed them in an open universe and how the in- 
falling particles could be given nonradial orbits. In the WZ 
model the initial specific binding energy of a shell is per- 
turbed in the following way: 



SE(M) = 5E, 



°\MJ 



(1) 



where M is the mass initially inside the shell. A halo with a 
fiat circular velocity curve near the centre results for e = 2/3. 
We use this particular case to generate our initial conditions. 

2.1 Collisionless component 

In the model of WZ the equations for the evolution of the 
collisionless component are: 



d 2 n 
dt 2 



(r, 2 +a 2 )^ 






t = GM * K *t 



GMiKi = 




t>U 



**=E 



(2) 
(3) 

(T) 
(5) 



where m;, Vi, J; and £; are the masses, radii, specific angular 
momenta, and turnround times of the shells, Mi is the total 
mass within shell i, and e, the eccentricity of the final orbit, 
is assumed to be the same for all shells (see WZ for more 
details). 

In order to avoid a divergence of the gravitational accel- 
eration near the origin the potential is softened in the usual 
way. By choosing an appropriate value of e we also give each 
shell enough angular momentum to prevent it from getting 
very close to the center. In most of our calculations we adopt 
e = 0.7 which results in almost isotropic velocity dispersions 
in the virialised regions of the halo. 

We integrate equations (2-5) by means of the force- 
average method. First of all we predict new positions and 
velocities using the accelerations calculated at the beginning 
of the time-step fi n : 



= n n + AoMt Vi n + 0.5 At DM fi 
— v t n + At D M fi n ■ 



(6) 

(7) 



Now, with the predicted positions we can calculate accel- 
erations at the end of the time step fi* and use a centred 
scheme to get positions and velocities for the beginning of 
the next step: 

Vl n+1 = Vi n + 0.5 Ainu (fi n + fi*) , (8) 

n n+1 = n n + 0.5 Atom (vi n + Vl n+1 ) . (9) 

Finally we just set fi n+1 = fi* and the step is complete. 
With this scheme we only have to calculate the forces once 
per timestep. The timestep AtoM is taken as: 



At i 



min Cd 



(10) 



where Cd m = 0.1, so that no individual particle moves more 
than 10% the softening length in a timestep. 

2.2 Gas component 

The initial conditions for the gas are generated in the same 
way as for the dark matter but with no angular momentum. 
We integrate the usual hydrodynamic equations, which in 
spherical symmetry can be written in the Lagrangian form, 
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dm 2 

— - = 4?rr p , 
or 

dv_, 2 dP 4tt d(r 3 Q) GM t r 

dt dm r dm • - 



du 
~dl 



-4tvP 



(r 2 +a 2 ) 2 



d(r 2 v) 3Q fdv d(r 2 v) 



(11) 
(12) 



dm 



■ + 



2 p \dr 



dr* y+-( r - A )>( 13 ) 



2.3 The cooling 

We consider a primordial plasma (i.e. hydrogen and he- 
lium only) with a helium mass fraction Y = 0.24. The 
integrations are started at redshift Zi — 500, just after 
recombination. The temperature of the gas is assumed to 
be uniform and equal to the temperature of the CMB, i.e. 
T = 2.726(1 + Zi) K ([Mather et at. 1994J). The initial frac- 










if 






<0 



otherwise 



P = ( 7 - l)pu ■ 



(13) 



(14) 



tional abundances relative to the total number density of 
hydrogen are taken to be: 

O 1 ' 2 

-5 "o 



{hili 



where p, v, P and u are the gas mass density, velocity, pres- 
sure and specific internal energy of the gas, A the radiative 
cooling rate, V the photoionization heating rate, 7 the ra- 
tio of specific heats, and Mt the total mass (gas and dark 



sp 
rti 



matter! inside the shell at, radius r. Tn order to treat shocks 



accurately we have implemented the Te nsor A rtificial Vis- 
cosity scheme of Tscharnuter & Winkler ( |l979[ ) (TAV). The 
artificial viscous pressure Q is only non-zero in regions that 
are undergoing compression (i.e. the divergence of the veloc- 
ity field is negative), and converts bulk kinetic energy into 
internal energy over the scale height I, which is taken to be 
I — XAr, where Ar is the radial thickness of the shell. The 
parameter A has to be big enough to damp out post-shock 
oscillations but small enough to have a good shock resolu- 
tion. We have found that a value A = \/6 works well and 
smears the shock over six zones. 

We set up a logarithmic grid in which equations (11- 



14) are discretized in the way prescribed by Benz ( 1991 ) 
(for details see Forcada-Miro 1997). Time integration is 
performed explicitly, except for the radiative processes, for 
which we switch between an explicit and an implicit method 



(see 2 3) Tn this way the timestep never becomes too smal 
when strong cooling occurs, but still, the limitations on the 
timestep for the gas are stronger than for the dark mat- 
ter. Consequently, we allow different timesteps for the two 
components. The timestep for the gas is taken as 



At g = min{AtcFL, At D M) 
where 

AtcFL = min Ccfl 



Ti —n-i 



i-i| + ^7(7 - l)ui 



(15) 



(16) 



and At dm is given in equation (10). In the simulations we 
use Ccfl — 0.25. The way we couple the two components 
is by advancing the dark matter particles a timestep AtDM , 
and then repeatedly advancing the gas until a further step 
Atg would take the gas ahead of the dark matter. At this 
point we just take the timestep required to synchronize the 
two components. The contribution of the dark matter to the 
total mass M t inside each gas shell is linearly interpolated in 
time. We use a force-average scheme similar to equations (6) 
to (9) to integrate the gas equations, but in this case the final 
advanced quantities are used to re-evaluate the accelerations 
at the end of each step before starting the next one. Forces 
are thus computed twice per time step. 
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(17) 

%e+ = %e++ = . (18) 

Equation (17) is the f it to t he residual ionization fraction of 
hydrogen by Peebles (1993). Helium would have rec ombined 
at an earlier epoch, so we can assume it neutral (Shapiro, 
Giroux fe Babul 199^ ). In all the simulations presented in 
this paper we took h — 0.5 and 0, o = 1. 

We use the ionization and recombination coefficients, 
and the cooling and heating rates revised by Abel (private 
comunication) . In order to reduce the computational time 
we interpolate the coefficients and rates out of a table con- 
structed beforehand. One has to bear in mind that there are 
many physical processes that can affect the emissivity of a 
hot primordial plasma in a drastic way (e.g. the presence of 
ionizing radiation, thermal conductivity by electrons, molec- 
ular hydrogen formation) and it is beyond the scope of this 
paper to study all their effects. 

We have assumed that electrons and protons have the 
same temperature. At the shock front it is the ions, and 
not the electrons, that are initially heated; the electrons are 
heated later by some mechanism which is at le ast as fast 
as Coulomb interactions with the positive ions (Shapiro & 



Moore 1976). The Coulomb equilibration time is given by 
Y 



te-iUH 



9.5 



1. 



1. - 0.75Y 



T 



10 5 K 



yr cm 



( [gpitzer 1962 ), and the cooling time, defined as 
_ {3/2)nkT vir 



COOL f 


^v,-' vir ) 


an be 


written 


tcoolTiH '' 


^ 6.6xio 4 
~ (i-y)M _ 




x r t 




'- L 10 B K 



more 

A(T) 



conveniently 

-1 



erg cm- J s- 



(19) 



(20) 



(22) 



yr cm 



where p is the mean molecular weight and Y is the primor- 
dial helium fraction by mass. In Figure la one can see that 
the Coulomb equilibration time is smaller than the cooling 
time up to very high temperatures, of the order of 10 8 K. It 
is thus a reasonable assumption to adopt a single post-shock 
temperature for electrons and ions. 

Whenever a mostly neutral plasma gets shock-heated 
up to the virial temperature T v i r , the ionization time, 
tion — (n e a)~ , can be longer than the cooling time {n e is 
the electron number density and a is the collisional ioniza- 
tion rate). The implication is that ionization of the plasma 
lags behind shock-heating, and line cooling can be enhanced 
by the higher neutral fractions. As one can see in Figure la, 
this is the case up to temperatures of the order of 10 6 K (i.e. 
V c ~ 160 km s _1 ). On the other hand, the UV flux radiated 
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by the hot shocked plasma can preionize hydrogen an d he- 
lium u pstream o f the shock provided V s > 110 km s _1 ( phull 
& McKee 197!:), where V s is the shock velocity in the rest 
frame of the shock. In our systems V s ~ 1.2 V c , so we expect 
the plasma to be preionized for halos with V c > 100 km s _1 . 
We do not take this effect into account in this paper. 

In higher circular velocity halos, the postshock heated 
gas quickly relaxes to the ionization equilibrium abun- 
dances. As the plasma reaches higher densities the cool- 
ing time decreases and eventually the gas begins to 
cool efficiently, while maintaining a higher ionized frac- 
tion. The reason is that the recombination time-scales, 
tree = (n e f3)~ , are much longer than the cooling time-scale 
(/? is the recombination rate). In Figure lb we can see that if 
the gas gets heated above 10 6 K it will cool out of ionization 
equilibrium. 

Because of the various possibilities, we do not assume 
equilibrium abundances for the different ions, but rather 
solve the full rate equations, 



1 


d %+ 


x e 


dt 


1 


d %c+ 


x e 


di 


1 


dl Hc++ 



1 s 



7j p g {/3h°£h° — a H+ a; H +} , 

£ Pg {/3ho°:ehc° + Qh c ++Ihc 
- (/3ho+ + "hc+ ) #He+ } 

1 a 



(23) 

(24) 



dt 



= £ P 3 {Ph C + X Hc+ -"Hc++%c++} > ( 25 ) 

where Xi is the fractional abundance of the species i relative 
to hydrogen, and j3i (eti) are the ionization (recombination) 
rates in units of cm 3 s _1 . The time t is given in the units, 

T = Ho' 1 = 9.78 h' 1 Gyr , (26) 

the density p~ g in the units, 

O TT 2 

v = : T^77 = 3 - 753 x 10" 2 V gr cm" 3 , (27) 

47rG 

and the factor C has the units of a reaction rate, 

C = 1.429 x 10 -13 (1 - Yp)" 1 h em's" 1 . (28) 

Equations (23)-(25) are written under the assumptions that 
the coronal limit is valid and the total number density of 
hydrogen varies on a much longer time scale than the ioniza- 
tion state of the plasma. It might happen that the problem 
is stiff, and then, cannot be handled efficiently by many in- 
tegration schemes. We use the subroutine LSODAR, coded 
by L. Petzold and A. Hindmarsh to dea l with stiff ordinary 
differential equations ( Hindmarsh 1983 ). 

If cooling/heating terms in the energy equation (equa- 
tion 13) are treated explicitly, the gas timestep has to be 
limited as well by the cooling time (equation 20). However, 
the timestep will become to small whenever strong cooling 
occurs. To prevent this fr om happening, th e code uses an 
operator splitting scheme ( Press et al. 1992 ). In a real case 
the microphysical processes operate continuously while the 
flow is evolving, but there is no way of maintaining this si- 
multaneity without paying a big price in the total computing 
time. The temperature at the end of a timestep T" +1 is the 
solution to the nonlinear algebraic equation, 



= 2(7-1) 

i 



ra+1,- »+l 



T" + 1 ) 



[« + 1 )ad + 



+1 QfX " +1 



T" + 1 ) Aig" 



(29) 



where (u™ + ) a d and p gi are the specific internal energy 
and the density which result from integrating the gas dy- 
namics equations without including the radiative terms, p 
is the mean molecular weight of a fully ionized primordial 
plasma, and T vlv is the virial temperature of the hosting 
halo, 



T vir = 6 x 10° Ho 
The factor L. 



Vc 



100 km s" 



K 



C = 2.377 x 10"^ (1 - Yp) 
and the function S(p g , T), 



Vc 



100 km s" 



S(p g ,T) 



r(pg,T)-A(p g ,T) 



n-a z 



(30) 



(31) 



(32) 



are both in units of erg s _1 cm 3 . 

The main complication to find a solution to equation 
(29) is that the ionisation state of the plasma depends on 
the final temperature. We have designed an iteration scheme 
to find the temperature and the ionisation state at the end 
of the timestep. A first guess for the temperature is taken 
to be TJj = T™. The rate equations are integrated, and 
the obtained fractional abundances are plugged in the right 
hand side of equation (29) to get a predicted temperature 
T* . If the error function, 

T* 

(33) 



J- 1 ) 



T, 



(i) 



is smaller than 0.001 in absolute value, then T™ + = T* and 
the integration is complete. If this is not the case a second 
guess for the temperature is taken to be, 



lf>= l-O.lsigCef 5 ; 



-COi 



-.(1) 



(34) 



where sig(e- ) is the signature of the error function. The rate 
equations are integrated and the whole process is repeated 
until the error function is smaller than 0.001, or changes 
sign. In the later case the root to equation (29) is bracketed 
with a 10% error, and the Va n Wijngaarden-Dekker-Brent 
root finder (Press et al. 1992) is used to find T™ +1 with a 



relative error smaller than 0.1%. 



2.4 Inner boundary 



Once the gas has been able to cool efficiently it piles up 
near the centre, and, since the calculations are performed 
in a Lagrangian frame, the nodes of the grid get extremely 
concentrated. This causes the code to require very small 



timesteps. Thoul & Weinberg (1995) got around this prob- 
lem by freezing at a certain distance r c those zones for which 
the cooling time has become much shorter than the dynam- 
ical time. However, they found that the amount of cooled 
mass was sensitive to this central boundary condition, and 
could vary by ~ 10% depending on the number of shells 
used. In addition, a recent numerical study of the stability 
of radiati ve shocks indicates that the se zones play an impor- 
tant role ( ptrickland fe Blondin 1995 ). Therefore we keep the 



innermost zone dynamically active, but we merge overlying 
zones with it whenever their temperature becomes smaller 
than a T c = 8000 K. 
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3 TESTS OF NUMERICS 

3.1 Adiabatic collapse 

Bertschinger (19851) found similarity solutions for secondary 
infall and accretion onto a seed mass in an Einstein de-Sitter 
universe. He studied not only collisonless and collisional col- 
lapses, but also the collapse of non-self-gravitating gas mov- 
ing in the potential of a collisionless component. These solu- 
tions provide a good testing ground for the dynamics of our 
code in the absence of the cooling processes. In WZ's model, 
Bertschinger's solutions correspond to e = 1 and no angular 
momentum (i. e. e — 1). However, in order to avoid inte- 
gration problems near the centre we prefer to use slightly 
non-radial orbits for the dark matter with an eccentricity 
e = 0.9. 

The self-similarity arises from the fact that there are 
no preferred scales in the initial perturbation, and therefore, 
the whole evolution of the system is determined only by the 
rate at which matter turns around. Following Bertschinger 
we nondimensionalize variables as follows: 



r t X 



v(r,t) 



T t 



V(\) 



p(r,t) = p H D(X) 
P{r,t) = p H (^) 2 p{\) 

m(r,t) = -np H r t 3 M{\) 



(35) 
(36) 

(37) 
(38) 

(39) 



where pu and rt are the mean cosmic density and the turn- 
round radius at time t. 

We first performed a calculation assuming a collision- 
less component and using 15000 shells. This simulation 
settled down to self-similar evolution which reproduced 
Bertschinger's solution almost perfectly except for small dif- 
ferences in the inner regions resulting from the fact that our 
shells have slightly nonradial orbits. The total energy was 
conserved up to 0.3%. We also calculated a collisionless col- 
lapse with e = 2/3 and e = 0.7. Our numerical s olutio n 
agrees well with that given by White & Zaritsky (1992). 
The mass profile at redshift z is linear in radius out to, 

Vc 



120 (l + z)" 3/2 h~ 



100 kn 



kpc 



and contains a mass, 



M vir ^2.8 x 10 11 (l + «) 



-3/2 



h~ 



V c 



(40) 



M e ■ (41) 



. 100 km B" 

The self-similar solution for a collisional gas is obtained 
using a grid with 500 particles. We have simulated only the 
adiabatic case 7 = 5/3, and find very good agreement with 
Bertschinger's similarity solution (see Table 1). The TAV 
scheme we have adopted proved to be very efficient in pro- 
ducing a sharp shock profile with practically no post-shock 
oscillations. We have also calculated a collisional collapse 
with both gas and a gravitationally dominant collisionless 
component (Qb — 0.05). In this case the agreement with 
Bertschinger's solution is also good (see Table 2), but we 
observe some post-shock oscillations, especially in the ve- 
locity profile. This is a spurious effect due to discrete jumps 
in the gravitational forces which occur whenever a dark mat- 
ter shell crosses a gas zone. To reduce such effects we have to 
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make sure that the shell masses are small enough compared 
to the gas zone masses. For equal masses this condition can 
be written as, 



Ndm ^ Qdm 



(42) 



gas zones respectively (Thoul & Weinberg 1995). The total 



energy was conserved up to 0.2%. 



3.2 Radiative processes 

In order to test our radiative processes subroutine we have 
calculated the thermal evolution of the intergalactic medium 
(IGM) in the absence of any per turba tion. We assume, fol- 
lowing Miralda-Escude & Rees ( 1994 ), that the IGM was 
fully ionized at a redshift z; n , and heated up to a tempera- 
ture Ti n . In their calculations these authors took Qb — 0.05, 
and we have also studied the cases Qt, = 0.01, 0.025 and 0.1. 
In the first case we find a good agreement with the pub- 
lished results. For a lower J7b the final temperature of the 
IGM turns out to be lower, while for a higher baryonic con- 
tent the final temperature is higher. This is the expected 
behaviour since the recombination rates decrease with den- 
sity, and consequently photoionization is not so effective at 
heating the gas. 

We have also studied the thermal evolution of an over- 
dense region, for which the density is given by the top-hat 
model until virialization, and is assumed constant there- 
after. Once again the agreement with Miralda-Escude & 
Rees (1994) is good. Now the final temperature of the cloud 
decreases as we increase the baryonic content of the Uni- 
verse. This is because at higher densities photoionization is 
not effective anymore and the dominant process is line cool- 
ing, which has an emissivity that increases with density. 



4 COLLAPSE SIMULATIONS 

4.1 Nonradiative collapse 

First of all we obtain the similarity solutions that arise from 
a power-law initial perturbation with index e = 2/3 in the 
absence of radiative cooling. If variables are scaled as in 
equations (35)- (39) the partial differential equations which 
describe the gas hydrodynamics transform into a set of or- 
dinary differential equations, 



(V - X)D' + DV' + -DV 
A 



2D = 



{V -\)V' = 



p' 



D 



2M_ 
9P 



(V-X)[j-^) =2(1-,) 
M' = 3A 2 £> , 



(43) 



(44) 



(45) 



(46) 



where primes denote derivatives with respect to A. 

Far away from the origin the solution for cold accretion 
applies, in which case analytic expressions can be obtained 
(for details see Forcada-Miro 1997), 



A = sin' 



G) 



9 — sind 



(47) 
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V = - cot - , 
2 \2 ' 



M 



3tt\ 2 TV 



(?) 



— sin# 



(48) 
(49) 

3V4 7 A(A-V)' (50) 

where 9 £ [0, 2tt]. Turnaround corresponds to a value of the 
parameter 9 t = n. 

At a given A = A s the nondimensional strong shock 
jump conditions, 



D = 



1/3tt\ 2 sin 2 (|) 



:(t) 



> - - As - ( — 4 ) (Vi - As) 



7+1 



Dz 



7 + 1 



7" 



Di 



7+1 
M 2 = Mi , 



£>i(Vi-A s ) 2 



(51) 
(52) 

(53) 
(54) 



provide the postshock values from which the ordinary differ- 
ential equa tions (43)- (46) are solved using a Runge-Kutta 
integrator (Press et al. f992). The solution obtained for 
7 = 5/3 is plotted in Figure 2. There is only one position 
of the shock A° fts 0.29 for which the boundary condition 
V(0) = M(0) = is satisfied. In this case the solution has 
an asymptotic behaviour near the origin, 



Z>~ A~ 



A" 



M~ A 



(55) 



The density profile is due to the fact that shells settle at a 
constant fraction of their turnaround radius. The pressure 
profile is a consequence of pressure forces balancing gravity 
in the static regions. If the shock is at a smaller A s then 
gravity rapidly dominates over pressure gradients and the 
fluid approaches adiabatic free-fall as it flows inwards: 

A" 1/2 , D~A- 3/2 , 



V 



D 1 



(56) 



Each fluid particle crosses the origin at a finite time, and 
thus, there is a mass flux through the inner boundary; in 
the limit A s = the solution for cold accretion is recovered. 
On the other hand, if A s > A° the mass goes to zero and the 
density diverges at a position < A D < A s . Notice that only 
the solution with A s = X° is physical. 

The introduction of a dominant collisionless component 
does not change the solution significantly from the purely 
gaseous case. We have performed the calculations for a bary- 
onic content J1b = 0.05 and eccentricities e = 0.9, 0.7, 0.5. 
In the self-gravitating case the shock position is only of the 
order of 1% smaller, and thus, the evolution in redshift of r s 
is well approximated for any baryonic content Ob G [0, 1] 
by the expression, 



r s ^I30(l + z)" 3/2 ft" 1 
and the gas mass within, 



M B 



3 x 10 u (l + z) 



-3/2 



100 km : 



Q B h' 



kpc 



Vc 



100 kn 



Mr. 



(57) 



(58) 



The mean overdensity enclosed by the shock is th en D s ~ 
120, s maller than the usually assumed value of 200 (White & 
Frenk 1991). It is also interesting that the infalling gas hits 



the shock with a velocity very close to the circular velocity 
of the inner halo. 



4.2 Radiative collapse 

The similarity solutions of the last section assume no radia- 
tive cooling. However, if the virial temperature of the system 
is above the sharp cut-off of the cooling curve at approxi- 
mately lO 4 -/!' the post-shock gas will be able to radiate, and 
the evolution of the system will be modified. According to 
equation (29) cooling should be important in halos with a 
circular velocity V c > 17 km s _1 . We have calculated the 
time evolution of systems with circular velocity 20, 40, 60, 
80, 100, 140, 180, 220 and 300 km s" 1 , and with gas frac- 
tions, Q.b — 0.05 and 0.1. Throughout this paper we adopt 
H = 50 kms" 1 Mpc" 1 . 

Radiative cooling introduces a new scale length into the 
problem, the cooling radius, defined as the radius at which 
the cooling time equals the lifetime of the system. The self- 
similarity we had in the nonradiative collapses is therefore 
broken. At a given epoch, the cooling radius and the shock 
radius are at different fractions of the turnaround radius 
for each circular velocity. For a given circular velocity these 
fractions change with time (see 4.2.1 and 4.2.2). 

The details of the behaviour of the radiative shock de- 
pend on whether the hosting halo has a circular velocity 
above or below 100 km s _1 . We therefore separate discus- 
sion of these two cases in the next two subsections. In the 
high circular velocity case the virial temperature is bigger 
than 3 x 10 5 K and thermal bremsstrahlung dominates the 
cooling. Line emission by H and He + provide the cooling 
processes in smaller halos. 



4-2.1 High circular velocity halos 

In Figure 3 we show the evolution of the radial structure of 
a halo like that of the Milky Way, i.e. V c = 220 km s" 1 . The 
variables in this figure have been scaled as in equations (34- 
38); O is the temperature scaled to the virial temperature of 
the system (equation 29). For this example we have chosen 
Q b = 0.05. 

At early times, the density profile inside the shock is 
well fitted by a power-law exponent —1.5, the expected so- 
l ution for a coolin g flow in an isothermal potential well 
( Bertschinger 1989| ). At later times, the profile has the ex- 
ponent of the nonradiative model, —2, from the shock ra- 
dius down to the cooling radius, and then at smaller radii it 
bends to an exponent of —1.5. The pressure profile is always 
slightly steeper than the density profile, because residual 
departures from hydrostatic equilibrium and the transition 
from the noncooling to the cooling regime all cause the tem- 
perature to increase inwards. At early times, the exponent 
for the pressure profile is roughly —1.8 over the regions we 
resolve, while at later times it is roughly —2.3 from the shock 
to the cooling radius, bending to —1.8 at smaller radii. The 
mass profiles follows in the outer regions the adiabatic col- 
lapse solution, and in the inner regions is dominated by the 
mass M coo i that has cooled below 8000 K. At early times 
Mcooi is close to the mass of the black hole in the case of 
cold accretion, and decreases as the system evolves in time. 

As gas passes through the cooling radius, radiative ef- 
fects become important; nevertheless, its temperature and 
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density increase continually until its cooling time becomes 
shorter than the local dynamical time and catastrophic cool- 
ing takes place. The reason is that as the gas radiates and 
flows inwards, compressional heating dominates over radia- 
tive cooling. At the point where catastrophic cooling takes 
place, the material gets very cold and dense and piles up 
near the center. We cut our plots at this point because our 
code cannot resolve the structure of the cold central region. 

Another consequence of this inflow is that the gas is 
more centrally concentrated than the dark matter, which 
in principle could be a possible way of explaning t he larger 
baryon fractions measured in rich galaxy clusters (Babul & 
Katz 1992). If we plot the baryon enhancement at a given 
time for different halos we see that the bigger the halo the 
closer we are to the adiabatic solution (Figure 4). On the 
other hand, for each individual halo the global expansion 
of the Universe reduces all densities and so makes radiative 
cooling progressively less efficient with time. As a result the 
gas enhancement gets weaker as the system evolves. This ex- 
plains why the profile of enclosed gas ma ss in Figure 3 gets 
steeper with time. However, White et al. (1992) proved that 
even assuming the most extreme case of cold accretion one 
cannot account for the observed baryon fractions in clus- 
ters of galaxies. Cooling processes alone cannot therefore 
enhance the cluster baryon fraction in an Einstein-de Sitter 
Universe to the values required for consistency with £1b as 
inferred from Big Bang nucleosynthesis. 

In Figure 5 we present the evolution of the shock radius 
and the cooling radius for circular velocities 140, 180, 220 
and 300 km s _1 , and in Figure 6 the evolution of the mass 
within these radii. In these simulations the central resolution 
is such that the shock forms at a redshift z ~ 50. At early 
times all the gas that goes through the shock is able to radi- 
ate away most of its internal energy. As a result, the shock 
wave and the cooling wave move together, their position in- 
creasing approximately as ~ t a , with a decreasing slightly 
with circular velocity. At later times the cooling wave sep- 
arates from the shock wave, because the post-shock density 
is no longer high enough for rapid cooling, and the shock 
position and the mass it contains grow in time in propor- 
tion to t 1 ' 2 , as in Bertschinger's (1986) similarity solution. 
The shock wave asymptotically approaches the nonradia- 
tive regime, its position and mass growing linearly in time. 
Figure 6 also shows the total mass that has cooled below 
8000 K, and so is included in the central zone. 



4-2.2 Low circular velocity halos 

The behaviour of the radiative shocks in systems with V c < 
100 km s~ x differs significantly from that in higher circular 
velocity halos. In order to clarify the situation in this case 
we have followed the evolution over longer timescales, cor- 
responding to epochs beyond z = 0. In order to do this it 
was necessary to reduce the central resolution slightly; the 
shock then forms at a redshift z ~ 25. 

The most noticeble feature of these simulations is that 
the shock front is unstable to large amplitude oscillations 
when 20 km s _1 < V c < 80 km s" 1 (Figure 7). The onset of 
these oscillations takes place earlier with decreasing circular 
velocity. Their amplitude and period grow as the system gets 
older, and, in the cases with 40 and 60 km s _1 , they die out 
at later times. The oscillations can also be seen in plots of 



the mass within the cooling and shock radii (Figure 8). It is 
interesting that during most of the oscillatory evolution no 
gas cools below 8000 K; only at the very end of the cycle 
does strong cooling allow a burst of mass accretion onto the 
central object. 

Radiative shocks have been found to be unstable to 
oscillations in other astrophysical contexts. Langer, Chan- 
mugam & Shaviv (1982) simulated time-dependent accre- 



tion onto white dwarfs and found that the cooling layer suf- 
fers oscillatory instabilities if the rate of accretion is above a 
certain critical value. Improved hydrodynamical simulations 
showed that these radiative shocks are unstable to oscilla- 
tions if the cooling rate A ~ p 2 T c has a dependence on 
temperature such as a < 0.6, in quite good agreement w ith 
linear stability analysis ([mamura, Wolf & Durisen 1984). 

These studies considered equilibrium ion abundances, 
and in our calculation we use time-dependent cooling. 
Therefore, our cooling function is dependent o n the his- 
tory of the plasma. Gaetz, Edgar & Chevalier (198C) ad- 
dressed the question of stability of radiative shocks with 
time-dependent cooling and concluded that, although a pre- 
cise stability limit cannot be established, simple power-law 
approximations of the cooling curve give a rough idea of the 
driving physics. 

With a power-law approximation the cooling function, 
the cooling time of a gas element goes as, 



tc 



(59) 



The post-shock temperature T s of the gas is related to the 
pre-shock infall velocity vi and the shock velocity v s through 
the strong shock jump conditions: 



(vi 



(60) 



If the shock moves outwards the gas is heated to a higher 
temperature, while the post-shock density drops (Figure 9). 
This leads to an increased cooling time, particularly for 
small c, and so to increased pressure in the region behind 
the shock. This excess pressure drives the shock further out- 
wards in a very short timescale. In this expansion phase the 
density, temperature and pressure profiles are rather fea- 
tureless, with slopes close to the cooling regime discussed in 
Section 4.2.1 (Figure 10). 

Eventually, the density of the downstream gas increases 
to the point where cooling again becomes substantial; pres- 
sure support for the shock then weakens and the shock slows 
down. The gas is now heated to a smaller temperature and 
can cool more efficiently (Figure 9) . As a result of this cool- 
ing instability a density peak forms behind the shock (Figure 
10). The rapid cooling of the gas in the density peak pro- 
duces a further loss of pressure support and the shock falls 
in with the flow. It is during this final collapse that cold gas 
is accreted onto the centre during a relatively short "burst" . 
A pressure build-up near the centre eventually causes the in- 
falling shock to strengthen again and to reverse its motion. 
As a result of the cosmic expansion, it takes longer for the 
post-shock gas to cool in each new cycle, causing the period 
and amplitude of the oscillations to increase. 

In our models we see such oscillatory behaviour when 
V c < 80 km s _ . This suggests that the critical slope of 
the cooling function for instability is c ~ 0. For the pri- 
mordial element abundances which we assume, c is reduced 
below the value c = 0.5 expected for bremsstrahlung as he- 
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Hum line cooling enhances the radiation at lower temper- 
atures. Including heavier elements in our cooling function 
would cause the transition to oscillatory behaviour to move 
to higher V c . 

Another interesting result is that in these low circu- 
lar velocity halos the gas component is more concentrated 
for a given fraction of the turnaround radius than in bigger 
systems (Figure 4). However, observations of dwarf galaxies 



seem t o indic ate that they are dark-matter dominated (Ko- 
rmendj 1988). Therefore, one needs some physical mecha- 



nism to dilute the baryon component in the smaller systems 
if one wants to reconcile this th eory o f galaxy formation 
with observations. Dekel & Silk ( 1986 ) roughly estimated 
that a starburst could cause loss of most of the gas in a halo 
with circular velocity below a critical value of the order of 
100km s _1 . Their proposed mechanism is supernova-driven 
winds, since gas absorbs radiation emitted by young stars 
very inefficiently. It is quite uncertain how much mass is 
lost by this mechanism which must be effective well before 
most of the gas turns into stars. 



5 THEORETICAL MODEL 

We desire to present physical explanations for the evolution 
in time of the shock and cooling wave. In order to accomplish 
this, we have use d as a starting point the ideas set out in 
White & Frenk ( [L99l| ). They pointed out that due to the 
general expansion of the Universe the shock will go through 
two very distinct epochs: 



(i) Infall-dominated phase: Initially the density is high 
enough that all the gas that gets shock-heated is able to 
cool (i.e. icooi < t, where t is the lifetime of the system). At 
this time the supply of cold gas is limited by the infall rate 
rather than by cooling. 

(ii) Cooling-dominated phase: After a certain time the 
density becomes too low and the post-shock gas cannot cool 
efficiently anymore (i.e. t coo i > i). At this stage the amount 
of cold gas is regulated by radiative losses. 

These two epochs have been clearly seen in the numerical 
simulations presented in Section 4, and we will study them 
now in more detail. 



5.1 Infall-dominated regime 

During the infall-dominated phase the gas is able to cool 
so efficiently that it rapidly loses pressure support. Conse- 
quently it flows towards the centre, most of the cooling layer 
being covered while the matter is still travelling at the post- 
shock velocity. The thickness of the cooling layer is then 
roughly, 



^^cool ^ ^s ^cool 



(61) 



As the gas cools it is compressed to high density. The whole 
process can be regarded as an isothermal shock, and thus 
the compression factor is, 



P2 ~ M B pi 

where A4 3 is the Mach number of the shock, 



(62) 



Ms 



p m p v B 
-yk B T c 



1/2 



>1 



(63) 



and T c is the gas temperature, which is of the order of the 
cut-off temperature in the cooling curve (i.e. T c ~ 10 4 K). 
The thickness of the cold layer can be roughly estimated 
from mass conservation, 



Ar co i d w M B 2/3 r s 



(64) 



Therefore, the shock will be at a position, 
r s « Ar C ooi > A?- C oid , (65) 

which implies that the shock moves in the quasi-equilibrium, 
tcooi ~ iff , (66) 

where ts is the freefall timescale, 



iff = 



1 / .W.. > ' J 

V24ivGp 



(tr<- 



-3/2 



(67) 



and p is the mean density of total mass within a given ra- 
dius. If the cooling time is bigger than the free-fall time 
adiabatic heating will increase the temperature of the gas 
thereby reducing the cooling time. In the opposite situa- 
tion the temperature rapidly decreases, and thus, the cooling 
time increases. 

At very early times Comp ton scatte ring provides the 
dominant source of energy loss ( pilk 1977 ). In this case, the 
cooling time depends only on the redshift, 



tco 



2.4 x 10 12 (1 + z) yr 



(68) 



According to equation (66), the evolution of the shock is 
qualitatively given by, 



-JT ~ ( 1 + z ) 



(69) 



The gas flows inwards until it becomes self-gravitating, 
which implies that M$ is roughly a constant, 



M s « D. B Mo 



(70) 



where M is the nondimensional mass ***what does this 
mean?*** of the black hole in the cold collapse solution. 
Therefore, the position of the shock evolves in time as, 

r s ~ Q^t 19 ' 9 . (71) 

The rate of energy loss due to Compton scattering de- 
creases with a high power of the redshift, and eventually 
two-body processes will dominate the radiative cooling of 
the system. In this case the cooling time depends on the 
nondimensional hydrodynamic variables as, 

tcooi ~ Di' 1 Vi 2 < 1 - e) (1 + z)- 3 , (72) 

and the qualitative evolution of the radiative shock is given 
by, 



*(A;c,n B ) = -DiVi 



-2(1- C )/Me\ 



-1/2 



(1 



-3/2 



(73) 



where the function ty(\;c, J7b) depends on two parameters: 
the slope of the cooling function c, and the baryonic content 
of the universe Qb- For a self-gravitating gas one has the 
asymptotic behaviour, 



*(A< l;c,fi B )~#o(fiB) A 1 



(74) 
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where ^(Ob) increases as the baryonic fraction decreases, 
and diverges for £3b —* 0. Therefore, the physical position of 
the radiative shock has a power-law behaviour in time, 



a (t;c, SIb) ~ r (Q B ) i 1 



(75) 



where r (flB) decreases with the baryonic fraction, and 
r (fi B ) -> as fi B -> 0. 

The radius of the shock wave in the infall-dominated 
regime grows faster than linearly in time; as a result the 
nondimensional shock position A s increases as the system 
evolves. The slope of the function *&(A; c, Qb) increases with 
A, and thus, the power-law behaviour of r s decreases ap- 
proaching a linear growth in time (Figure 11). Eventually, 
the postshock cooling time becomes bigger than the lifetime 
of the system, and the shock enters the cooling-dominated 
regime. An upper limit to the position at which this tran- 
sition happens is provided by the nonradiative shock (i.e. 
A s ~ 0.29). At this point, r s has roughly a power-law be- 
haviour, 

1.5 -c 



(76) 

Therefore, if one fits a single power-law to the time evolution 
of the radiative shock one should obtain a value in between 
these two limits. 

One has to keep in mind that the softened gravity pre- 
vents the collapse beyond the resolution limit of the code, 
and the shock forms at a position close to the softening pa- 
rameter a. For r < a a cold dense core forms, where the den- 
sity is constant in radius and the material is at rest. At this 
stage the shock strongly violates the quasi-equilibrium con- 
dition, with £ coo i <C iff. In the numerical simulations which 
do not present an oscillatory instability we have fitted an 
analytical expression of the form r s = a + bt a to the evolu- 
tion in time of the shock radius (Figure 5 and Figure 7). The 
fitted exponent is in all cases between the predicted limits, 
and for the high velocity halos is fairly constant, with a value 
close to a — 2. This is the exponent expected when Compton 
cooling is still a significant source of cooling (Figure 12). 



5.2 Cooling-dominated regime 

The transition to the cooling-dominated regime happens 
when the postshock cooling time equals the lifetime of the 
system. According to equations (66) and (69) the mean total 
overdensity has a well defined value at this time. This is in 
agreement with the fact that in the simulations the nondi- 
mensional position at which the shock transits is the same 
for all halos with a very small scatter (Table 3). Also the 
postshock temperature is very close to the virial tempera- 
ture of the system. 

The nondimensional position of the shock at the time of 
transition is very close to the case 7 = 4/3. In fact, radiative 
cooling weakens the pressure support of the gas against grav- 
ity, making the equations of state less stiff; in other words, 
radiative cooling leads to an effective ratio of specific heats 
7 c ff < 5/3, and, if the gas cools strongly, even 7 B g < 4/3. 
The numerical simulations suggest that one can think of the 
evolution in time of a collapsing halo as "passing" through 
a series of adiabats. At early times, when radiative cooling 
is very strong, one has 7 e ff < 4/3. This exponent increases 
with time towards 7 = 5/3. The value 7 = 4/3 delimits 
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the boundary between the infall-dominated and the cooling- 
dominated regime. 

With the appropriate choice of temperature and gas 
density one can obtain a fairly good estimate of the redshift 
Zd at which the cooling-dominated phase starts. If one takes 
a temperature equal to the virial temperature of the halo, 
and a density equal to the mean gas overdensity inside the 
nonradiative shock the predictions compare well with the 
results from the simulations (Figure 13). Now, for a given 
halo, the physical size at the time the transition takes place 
is simply, 

-1 -3/2 
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10 



1 



' 2rf 



1 



ft" 



Vc 



100 km s~ 



kpc 



(77) 



where z* { = 3 is a characteristic redshift, and the total mass 
of cold gas enclosed in solar masses is, 
-1 -3/2 
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n B h~ 



Vc 



100 km s" 



(78) 



The simulations show that as soon as the shock enters 
the cooling-dominated regime a hot halo starts building with 
an isothermal density profile. The temperature has a slight 
increase inwards (Figure 3), but we will assume it is a con- 
stant and equal to the virial temperature of the hosting halo. 
At any given time there is a radius r c < r s at which the gas 
density is high enough so the cooling time equals the life 
time of the system, 

Pg(^s) 



Pg(Tc) 
and since the halo is isothermal, 

1/2 
r c (t) « r s (t) 
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£cool(f sj 



(79) 



(80) 



Once the shock has entered the cooling-dominated regime, it 
gradually approaches the nonradiative behaviour and grows 
linearly in time. At this stage the postshock cooling time 
grows proportional to t 2 , and then, 

r t i 1/2 

r c (t)«r.(to) =— , (81) 

_ cool _ 

where t is a reference epoch at which the shock presents the 
nonradiative behaviour, and t* ool is the post-shock cooling 
time at the reference epoch. 

The cooling time t* ool is related to the estimated time 
at which the shock enters the cooling dominated regime £ c f 
as. 



tc 



El 



(82) 



where D s is the mean gas overdensity inside the nonradiative 
shock and D s ~ 27 is the self-similar postshock density. 
Combining equations (81) and (82) we have obtained the 
following analytic approximation for the evolution in time 
of the position of the cooling wave, 
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(83) 



In Figure 5 and Figure 7 one can see that the analytical 
expression given in equation (83) is a rather good approx- 
imation to the position of the cooling wave once the shock 
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wave presents the nonradiative behaviour. Until this hap- 
pens the postshock density decreases in time faster than in 
the self-similar nonradiative evolution, so the cooling radius 
will be at a smaller fraction of the shock radius than pre- 
dicted by a power-law scaling. In fact, in Figure 5 one can 
see that in the early stages of the cooling flow the radius r c 
decreases with time. On the other hand, the cooling radius 
cannot move inwards faster than the flow, so it might hap- 
pen that the cooling wave is stationary at the early stages, 
with no mass accretion. The numerical code only gives the 
position of the cooling radius whenever a cell can first cool 
efficiently, and thus, a stationary cooling wave will produce 
a "gap" in the data. This "gap" is clearly seen for V c — 300 
km s _1 in both the radius (Figure 5) and the mass evolution 
(Figure 6), and barely for V c — 220 km s~ ■ This effect is 
stronger with increasing circular velocity because the power- 
law a is bigger (Figure 12). 

In the self-similar solutions of Bertschinger (1989) the 
mass accretion rate at the cooling radius is approximately, 



M c rs 4 7r r c p g (r c 



(84) 



so we model the rate at which cold gas accumulates by this 
simple inflow equation. Since the halo is isothermal down to 
the cooling radius, 



2 

r c ft 



(r c 



Pe(r s ) 



(85) 



and at later times, when the system is well into the cooling- 
dominated regime, the total mass accreted by the shock is, 

3 



M s « 4 ty r s d p e (r s ) . 

Now we can write equation (84) as, 



(86) 



(87) 



physical behaviour of the shocks, we have been able to build 
a simple model which explains very well the kinematics both 
of the shock wave and of the cooling wave in all cases which 
are not unstable. At early times the shock front evolves in 
such a way that the cooling time is of the order of the free- 
fall time. In the bigger halos Compton cooling dominates, 
and the radius of the shock grows proportional to t 2 . In 
smaller halos the slope is a function of the slope of the cool- 
ing function at the virial temperature of the system. The 
transition to the cooling-dominated regime happens at the 
same nondimensional position in all halos. This is close to 
the position of a nonradiative shock with a ratio of spe- 
cific heats 7 = 4/3. The time at which each system becomes 
nonradiative is well estimated assuming a temperature equal 
to the virial temperature of the system, and an overdensity 
equal to the mean overdensity inside the nonradiative shock. 
Whenever oscillations appear, our model provides only 
a rough guide to the behaviour. However, one should keep in 
mind that there are important processes that we have not 
taken into account, and that affect systems with circular 
velocities precisely in the range where the shock wave os- 
cillates. The cold dense gas that accumulates in the central 
regions will form stars, which will provide a flux of ioniz- 
ing radiation. In addition we should consider a metagalactic 
background of UV radiation due to high redshift quasars. 
The effect of such radiation is to maintain a higher ioniza- 
tion fraction which reduces the line cooling. We may expect 
the oscillatory behaviour to be absent in this case, and lower 
circular velocity dwarf galaxies be prevented from forming 
(Efstathiou 1992). We will consider such effects further in a 
future paper. They have convinced us that further analysis 
of the oscillatory solutions is unlikely to be worthwhile. 



According to equation (87), and since the shock evolves like 
in the nonradiative case, the mass inside the cooling radius 
grows in time proportional to t 1 ' 2 . In Figure 6 and Figure 
8 one can see that a simple power-law scaling, 



M c (t) w M s (t c f) 



1/2 



is a good approximation for the mass within the cooling 
radius. 
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6 CONCLUSIONS 

1-D simulations like the ones we have done are a powerful 
tool to enhance our understanding of the physics involved 
in the process of galaxy formation. In particular, we have 
concentrated on the shocks that appear in the collapse of 
a protogalaxy, and they have been proved to present the 
features already known from previous studies of strongly ra- 
diative shocks in other contexts. The shocks are subject to 
an oscillatory instability whenever line cooling dominates 
the radiative losses behind the shock front. Due to the ex- 
panding background in which the system is embedded, the 
period and the amplitude of the oscillations increase with 
time. 

We believe that we now fully understand the effects 
of radiative cooling in the shocks that appear in spherical 
models for galaxy formation. Despite the rather complicated 
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Figure 1: (a) Coulomb equilibration time (straight solid line), cooling time (solid), and collisional ionisation times for HI (dotted), Hel 

(short-dashed) and Hell (long-dashed) for a primordial plasma with Qg = 0.05, h = 0.5 and residual ionised fraction taken from 

Peebles (1993). (b) Cooling time (solid) and recombination times for HII (dotted), Hell (short-dashed) and Helll (long-dashed) for a 

primordial plasma with Hb = 0.05, h = 0.5 and equilibrium ionised fraction. 
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Figure 2: Self-similar velocity (top-left), density (top-right), pressure (bottom-left), mass (bottom-right) profiles for shocked accretion 

in the case e = 2/3 and 7 = 5/3. Solid line: solution satisfying the inner boundary V(0) = M(0) =0. Dashed lines: solutions with A s < A° 

or A s > A°, where \° is the position of the shock for the solution plotted in solid line. The expected asymptotic slopes are indicated. 
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Figure 3: Radial profiles at different redshifts for a halo with V c = 220 km s _1 . Top: Density. The two slopes plotted are -2 and -1.5. 

Middle: Temperature. The slope plotted is -0.3. Bottom: Mass. The top horizontal dashed line is the mass of the black hole in the case 

of cold accretion and the bottom one is the total mass that has cooled below 8000 K. The solid line is the self-similar mass profile in 

the case of adiabatic collapse. In all cases the left vertical dashed line is the position of the cooling wave and the right one is the 

position of the shock wave. 
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Figure 4: Mass of baryons to total mass ratio versus fraction of the turnaround radius at z = 0. In short-dashed line, from left to right 
one has halos with circular velocity 300, 140, 80, 40 and 20 km s~ 1 . In solid line one has the adiabatic collapse (bottom) and collapse 

onto a black hole (top). In all cases f2e = 0.1. 
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Figure 5: Evolution in time of the shock wave position (solid line) and the cooling wave position (crosses) as given by the numerical 

simulations for halos with V c > 100 km s _1 and Q^ = 0.05. Long-dashed lines correspond to the adiabatic evolution of the shock. 

Dotted-lines correspond to a growth in time proportional to t a , and an off-set equal to the softening parameter. The value of the 

exponent a is given in each panel. Short-dashed lines correspond to the analytical expression given in equation (6.31). 
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Figure 6: Evolution in time of the gas mass within the shock wave (solid line), the cooling wave (crosses), and the mass that has 

cooled below 8000 K (dotted line) as given by the numerical simulations for halos with V c > 100 km s —1 and Qg = 0.05. Long-dashed 

lines correspond to the adiabatic evolution of the shock. Short-dashed lines indicate a growth in time proportional to t ' . 
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Figure 7: Evolution in time of the shock wave position (solid line) and the cooling wave position (crosses) as given by the numerical 

simulations for halos with V c < 100 km s — 1 and Hb = 0.05. Long-dashed lines correspond to the adiabatic evolution of the shock. In 

the bottom right panel the short-dashed line corresponds to the analytical expression given in equation (6.31). The dotted-lines 

indicate a power-law t , with an off-set equal to the softening parameter. 
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Figure 8: Evolution in time of the gas mass within the shock wave (solid line), the cooling wave (crosses), and the mass that has 

cooled below 8000 K (dotted line) as given by the numerical simulations for halos with V c < 100 km s — 1 and Cb = 0.05. Long-dashed 

lines correspond to the adiabatic evolution of the shock. In the bottom right panel the short-dashed line corresponds to a growth in 

time proportional to 4 1 / 2 . 
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Figure 9: Evolution in time of the post-shock density (top) and the post-shock temperature (bottom) for a halo with V c 
and f^B = 0.05. Short-dashed lines indicate the times at which radial profiles are plotted in Figure 6.9. 
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Figure 10: Radial profiles at different redshifts for a halo with V c = 40 km s _1 and Ob = 0.05. Top: Density. Middle: Temperature. 

Bottom: Pressure. The left panel is a snapshot taken just after the shock reached maximum expansion. The right panel was taken 

during the expansion phase. The short-dashed lines indicate the power-law behaviour —1.5, —0.3, —1.8 for the density, temperature and 

pressure respectively. 
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Figure 11: Variation of the slope of r B (t) with A for different slopes of the cooling function c in the case Ob = 1- The long-dashed lines 

indicate the position of the nonradiative shock and the turnaround radius. 
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Figure 12: Solid lines: Predicted limits for the power-law evolution in time of the position of a radiative shock r s ~ t a in the case 

two-body process dominate the radiative cooling. The lower limit is the analytical approximation to the numerical limit (dotted line). 

The short-dashed line is the expected exponent if Compton cooling is the dominant cooling. Filled squares: Fitted exponent to 

simulations with Q^ = 0.05. Empty squares: Fitted exponent to simulations with Qb = 0.1. 
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Figure 13: Predicted redshift for the onset of a cooling flow at the central regions of a collapsing halo in a flat unvierse with a 
baryonic fraction Cb = 0.05 (solid line) and Q^ = 0.1 (dashed line). Filled squares: numerical simulations with f^B = 0.05. Crosses: 

numerical simulations with Hb =0.1. 
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